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Section  I 
PLASMA  THLORY 


There  is  no  report  for  this  Section  this  quarter.  The  remain- 
der of  Nevins1  work,  reported  last  quarter,  was  completed. 

Section  II 
SIMULATION 

A.  OH AS  I L INIAK  SIMULA! (ON  01  A MIRROR  MACHINE 
Dr.  Y.  Matsuda  with  Dr.  H.  L.  Berk  ( LLL ) 

While  work  continues  in  this  area,  there  is  no  special  report 
to  make  this  quarter. 

B.  SELF-HEATING  OF  Id  THERMAL  PLASMA 
A.  Peiravi  (Prof.  C.  K.  Birdsall) 

We  have  already  presented  simulation  results  for  the  self- 
heating times  of  a one-dimensional  electrostatic  thermal  plasma  and  a 
comparison  of  cost  of  simulation  vt>  heating  time.  We  now  present  com- 
plete results  of  self-heating  times  vM  wpeAt  and  Ap/Ax  ^or  nearest-  9ricJ 
point  weighting,  NGP,  linear  weighting,  CIC,  and  quadratic  spline  weight- 
ing, QS,  along  with  gain  considerations. 

The  heating  time  was  found  to  be  proportional  to  N^.  + Np  for 
a fixed  ratio  of  Ap/Ax  . Figures  1,2,  and  3 show  the  self-heating  times 

divided  by  N.  + N„  A.  / Ax  for  different  values  of  id  At  . The 

C D D pe 

dashed  line  drawn  through  the  figures  is  v (At/Ax)  ^ 3/2  for  NGP  and 
v^(At/Ax)  = l/2for  CIC  and  QS.  These  lines  indicate  the  optimal  choices 
of  parameters.  Figures  4,  5.  and  6 show  ratios  of  self-heating  times  of 
the  different  weights  vM  Ap/Ax  for  iopeAt  = 0.1,  0.2,  0.3  , respectively. 
These  figures  indicate  that  CIC  is  as  much  as  70  times  better  than  NGP, 
and  QS  is  as  much  as  650  times  better  than  NGP.  Actual  measurements  of 
cost  of  simulation  per  particle  per  time  step  on  the  CDC-7600  MFE  compu- 
ter at  LLL  show  that 
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TNGP  = 9 usec/par t i cle/t ime  step 

T =116"  " " 

CIC  ° 

T = oL  " ii  ii 

QS 


Defining  gain  of  using  a higher  order  weighting  scheme  as 

. increase  in  self-heating  time 

gain  = : r-2 

increase  in  computer  time 

and  going  through  the  optimal  path,  we  have  gains  as  in  Table  1. 


„ w At 
\ pe 

gain\ 

0.1 

0.2 

0.3 

CIC/NGP 

TJ-  20  8 

QS/NGP 

irf  = 72-9 

fl  S 

oo|o 

II 

v-o 

x- 

= 91’6 

QS/CIC  = 

( QS/NGP ) / ( C 1 C/NGP) 

6.1 

6.5 

3-0 

T.D1  c , Gain  in  going  to  a higher  order  algorithm.  The  ratios 

IMULt  I.  . .«  • . i/-.  . • /. 

inside  the  table  are  increase  in  self-heating  times  (i.e., 
reduction  in  error  in  energy)  over  increase  in  cost  (determined  on 
the  CDC  7600  MFE  computer  at  LLL). 


The  self-heating  time  can  be  increased  considerably  by  trunca- 
tion in  k-space.  Figure  7 shows  the  self-heating  time 

k /k.  . for  the  different  schemes  used  (keepinq  everything  fixed  but 

max  last  a ' 3 

varying  k]ast  for  each  scheme).  The  gain  in  self-heating  time  due  to 

k-space  truncation  is  almost  proportional  to  k /k,  for  NGP  and  is 

max  last 

close  to  but  not  quite  proportional  to  (k  /k.  )2  for  CIC  and 

max  laxt 

(k  /k,  )3  for  QS.  Thus,  k-space  truncation  further  increases  the  gain 

max  last  K 

of  CIC/NGP  and  QS/NGP  and  QS/CIC.  Table  2 shows  approximate  gains  with 
k-space  truncation. 
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v a)  At 
\ pe 

gain 
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*(k  /k.  ) 
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H.9 

20.8 

30. A 

QS/NGP 

x(k  /k  )2 

max  last 

72.9 

135- A 

91.6 

QS/CIC 

*(k  /k  ) 
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6.  1 

6.5 

3.0 

Gain  wi th 

TABLE  2.  k-space 

t runcat i on . 


C.  VELOCITY-SPACE  RING-PLASMA  INSTABILITY,  MAGNETIZED 
J.  K.  Lee  (Prof.  C.  K.  Birdsall) 

We  have  summarized  for  journal  publication  our  theory  and  simu- 
lations (two  papers)  on  the  linear  and  nonlinear  effects  of  this  short- 
time  scale  flute-like  instability.  This  instability  is  unrelated  to  colli 
sions,  resulting  from  the  interactions  of  energetic  ion  components  (ring- 
like in  velocity-space,  but  not  necessarily  monoenerget i c)  with  cooler 
Maxwellian  target  ions  and  electrons.  We  have  included  a wide  range  of 
electromagnetic  effects.  A third  report  on  the  1 } dimensional  hybrid 
code  (two  versions:  electrostatic  and  electromagnetic)  will  be  written 
separately  for  submission  to  the  Journal  of  Computational  Physics.  The 
table  of  contents  of  the  first  two  papers,  including  some  new  results 
obtained  during  this  quarter,  is  as  follows 

. 


VELOCITY  - SPACE  RING-PLASMA  INSTABILITY,  MAGNETIZED 

PART  I : THEORY 


I.  INTRODUCTION 

II.  ELECTROSTATIC  DISPERSION  RELATION 

III.  CLASS!  FI  CATION  OF  REGIMES  biy  RING  VENS1TV 

(/)  Weak  Ring  Regime  (5  x 10_t+  < R ' 5 x 10-3) 

( 2 ) Intermediate  Ring  Regime  (5xiO~,sRs5xIO~7) 

(3)  Strong  Ring  Regime  (5X10~7<R<5.0) 

IV.  EFFECTS  o{s  VARIOUS  PARAMETERS 

(/)  Add  ing  Electron  Dynamics 

(2)  Stabilizing  Effects  of  the  Plasma  Thermal  Spread 

(3)  Stabilizing  Effects  of  the  Ring  Thermal  Spread 

(4)  Effects  of  the  Target  Plasma  Density 

(5)  Effects  of  Unlike  Ion  Species 

V.  ELECTROMAGNETIC  MODIFICATION 

VI.  CONCLUSIONS 
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Theoretical  Estimations 
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D.  ONE-DIMENSIONAL  ESI  CODE  FOR  SHEATHS 
Yu-Jiuan  Chen  (Prof.  C.  K.  Birdsall) 

There  is  no  new  work  to  report  this  quarter. 

f.  PARTICLE  TRAJECTORIES  IN  A CUSP  FIELD 
Yu-Jiuan  Chen  (Prof.  C.  K.  Birdsall) 

Previous  investigation  by  Levine  e .t  cut.*  of  particle  loss  in  a 
toroidally  symmetric  cusp  shows  a reduction  in  loss  rate  due  to  symmetric 
geometric  effects;  this  is  because  the  conservation  of  canonical  momentum 
along  the  symmetric  axis  prevents  some  particles  from  penetrating  (and 
being  lost  in)  the  mirror-like  magnetic  field  at  the  cusp  point.  It  was 
pointed  out  by  Boozer  and  Levine2  that  in  a two-dimensional  cusp  (Fig.  1), 
a large  fraction  of  the  particles  in  the  field-plasma  boundary  move  in 
definite  guiding-center  orbits.  Therefore,  such  particles  will  remain 
trapped  between  Coulomb  collisions  in  the  absence  of  microinstabilities. 
The  trapping  of  particles  on  the  open-cusp  field  lines  is  similar  to  the 
particle  trapping  in  a mirror  machine. 

We  have  re-examined  the  particle  trajectories  in  an  axially 
symmetric  2d  cusp  by  using  Levine's  code  which  integrates  Newton's  equa- 
tion of  motion  using  the  fourth-order  Runge-Kutta  method  to  determine  the 
particle  motion  in  the  magnetic  field.  The  magnetic  field  is  given  analy- 
tically as  in  Ref.  2.  For  a sharp  boundary,  the  form  is  given  as 

(a)  if  x2/3  + y2/3  < 1 , Az  = Bx  = By  = 0 (I) 

(b)  if  x2/3  + y2/3>  1 , z = x + iy 

« 

'Morton  A.  Levine,  Allen  H.  Boozer,  G.  Kalman  and  Pradip  Bakshi,  "Particle 
Loss  in  a Toroidally  Symmetric  Cusp",  Phys.  Rev.  Lett.  28,  1323  (1972). 

2Allen  H.  Boozer  and  Morton  A.  Levine,  "Particle  Trapping  in  Magnetic  Line 
Cusps",  Phys.  Rev.  Lett.  3J.,  1287  (1973). 
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FIG.  1 The  configuration  of  the  2d  cusp.  §c  is  the  cusp  field  due  to 
four  line  current  wires.  §z  is  the  field  perpendicular  to  the 

cusp  field. 


FIG.  2 The  coordinates  v and  <)>  on  the  sphere  of  possible  velocity 


di rect ions . 
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? lx  I * ' J 

U = 2[z  + (z  - 1)J] 

L 2 H(U  + 4/U)  - 4z(U  - 4/U)_i  J 

B = 2/[(U  + 4/U)*  - (U  + 4/U-4L)*J 

Bx  = Re (B) 

By  = - I m (B) 

A = Im  | 3 (B2  + B'2)/81  (2) 

For  a sheath  with  a finite  thickness  A , the  vector  potential  near  the 
plasma  is  given  as 

Az(  *,y  ) = Bc  [ A - Atanh (A/ A ) ] (3) 


The  coordinate  system  used  for  describing  particles  is  given 
in  Fig.  2.  Figure  3(a,b,c,d,e)  show  the  trajectories  of  particles  with 
various  initial  conditions.  In  these  calculations,  time  is  measured 
in  inverse  cyclotron-frequency  units.  The  sheath  thickness,  A .was 
chosen  to  be  2.6  times  the  ion-cyclotron  radius,  p.  , where  p.  was 
choosen  to  be  0.02  times  as  large  as  the  distance  between  the  center  of 
the  plasma  and  the  cusp.  The  normalized  uniform  field  B is  1 and 
[Bz  + BZ]J  = 1 . The  equation  of  the  field-plasma  boundary  is 


X2/3  + y2/3 


(4) 


In  Fig.  3(a)  the  guiding  center  of  a particle  initially  on 
the  boundary  at  the  point  (0.356,  0.356)  with  v^/v  = 0.h  and  <f>  = 30° 
does  not  move  during  the  time  interval  (0.050,  826.437).  In  Fig.  3(b) 
the  particle  initially  is  at  2 cyclotron  radii  in  the  sheath  with  v^/v = 
0.6  and  $ = 45°  ; the  guiding  center  moves  down  one  side  then  returns. 

In  Fig.  3(c)  the  particle  starts  from  the  cusp  point  (1.0,  1.0*  10"(>)  with 
v^/v  = 0.2  and  <j)  = 45°  ; it  moves  up,  then  circles  four  sides.  Fig.  3(d) 
shows  another  type  of  trajectory  in  which  the  particle  moves  up  and  only 
circles  one  side.  In  Fig.  3(e)  the  particle  starts  from  the  position 
(0.396,  0.396),  i.e.,  3 cyclotron  radii  from  the  boundary  or  0.4  cyclotron 
radii  outside  the  sheath,  and  escapes  from  the  cusp  point  very  quickly. 


lb 


These  figures  show  that  particles  in  the  sheath  move  in 
definite  guiding  center  orbits  even  though  the  magnetic  field  has  large 
variations  in  magnitude  and  direction  in  the  sheath.  Because  the 
conservation  of  the  canonical  momentum, 


eA 

p = mv  + — - , (5) 

z z c 


due  to  the  geometric  symmetry  and  the  small  spatial  variation  of  Az  in 
the  sheath,  the  trapped  particles  always  remain  within  a few  cyclotron 
radi i of  the  plasma. 

In  Fig.  3(c)  the  axis  of  the  trajectory  is  shifted  slightly 
to  the  right  in  the  second  quadrant.  This  shift  raises  questions  about 
the  possible  jumps  in  the  usual  adiabatic  invariants 


J 


p dx 

1 


(6) 


(where  p^  and  x^  are  the  components  of  momentum  and  position  vectors 
perpendicular  to  the  total  magnetic  field),  the  magnetic  moment  p , and 
the  conserved  canonical  momentum  . These  questions  will  be  addressed 
in  the  next  quarter.  Also,  the  code  TIBRO-X  will  be  used  to  study  the 
particle  orbits  in  either  the  analytic  magnetic  field  or  the  self-consis- 
tent field  with  the  toroidal  effect. 


F.  PLASMAS  WITH  FIELD  REVERSAL 

Douglas  Harned  (Prof.  C.  K.  Birdsall) 

At  the  present  time  there  is  a great  deal  of  interest  in  field 
reversed  configurations  for  plasma  confinement,  partly  because  of  their 
attractiveness  for  potential  fusion  reactor  designs  as  well  as  for  under- 
standing space  plasmas  (e.g.,  the  bow  shock).  Before  beginning  work  in 
the  area  of  field  reversal,  it  has  been  necessary  to  determine  the  status 
and  direction  of  current  research  in  this  area.  Of  particular  laboratory 
interest  are  field  reversed  mirrors  and  intense  ion  rings. 
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Research  involving  the  field  reversed  mirror  has  been  aimed 
at  calculating  electron  return  currents  and  the  effects  of  quadrupole 
fields.  The  work  with  quadrupole  fields  has  been  approached  largely  by 
drawing  analogies  with  the  RECE-Berta  electron  ring  experiments  at 
Cornell.  Analytic  and  simulation  studies  are  in  progress  at  Livermore 
on  enhanced  resonant  diffusion  due  to  the  application  of  a quadrupole 
field.  These  results  have  so  far  given  only  very  rough  agreement  with 
RECE-Berta.  The  study  of  quadrupole  fields  is  now  being  applied  to  the 
field  reversed  mirror  by  their  inclusion  in  simulations  using  both 
TIBRO  and  SUPERLAYER  for  single  particle  trajectories. 

With  regard  to  electron  return  currents,  it  is  still  uncertain 
as  to  whether  or  not  these  currents  can  effectively  prevent  the  achieve- 
ment of  field  reversal.  Work  by  Baldwin  and  Fowler1  demonstrates  the 
need  for  some  process  (e.g.,  Ohkawa  currents)  to  prevent  neutralization 
of  the  current  at  the  field  null.  SUPERLAYER,  which  heretofore  has  ignored 
electron  effects,  is  now  being  modified  to  enable  the  investigation  of 
electron  return  currents. 

In  other  areas,  concerning  FRM,  stability  analyses  have  shown 
that,  in  order  to  avoid  MHO  and  tearing  modes,  the  axial  length  of  the 
plasma  as  well  as  the  major  and  minor  radii  of  the  field  reversed  region 
must  be  on  the  order  of  a few  Larmor  radii.  In  the  area  of  transport, 
very  little,  if  any,  work  has  been  done  relating  to  field  reversed  plasmas. 

At  Cornell,  ion  rings  have  been  investigated  (through  both 
simulation  and  theory)  from  the  standpoints  of  both  macro-  and  micro-insta- 
bilities. These  studies  have  determined  that  the  most  stable  ring  is  an 
intermediate  configuration  somewhere  between  a long  layer  and  a bicycle 
tire.  A number  of  codes  have  been  developed  to  study  ion  ring  behavior, 
the  most  recent  of  which  is  a 3-dimensional  hybrid  code  using  ring-like 
supe rpar tic les . 2 

'D.  E.  Baldwin,  T.  K.  Fowler,  "A  Reassessment  of  the  Requirements  to  Obtain 
Field  Reversal  in  Mirror  Machines",  UC1D-17691,  Dec.  20,  1977- 

2 A . Friedman,  R.  N.  Sudan,  J.  Denavit,  "A  Linearized  3'D  Hybrid  Code  for 
Stability  Studies  of  Field  Reversed  Ion  Rings  and  Mirror  Confinement", 

Bull.  Amer.  Phys.  Soc.  22,  9,  P-  1070,  1977. 
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Among  the  potential  areas  for  our  research  are  the  following: 

1)  Utilizing  TIBRO-X  to  investigate  particle  trajectories  in  field 
reversed  configurations.  Here  the  interest  could  focus  on  two  areas: 

a)  Trajectory  behavior  near  the  reconnection  points.  It  is  of 
interest  to  see  what  effect  the  reconnection  points  will  have 
on  confinement.  Although  this  region  of  the  plasma  lends 
itself  well  to  simple  models  (such  as  a reversed  current  loop 
about  the  center  of  a mirror)  it  may  prove  difficult  to  simu- 
late because  finite  time  steps  will  be  inaccurate  near  the 
reconnection  points. 

b)  Confinement  of  particles  in  closed  field  lines.  This  simulation 
requires  a much  more  complex  model  since  the  very  nature  of 

the  closed  region  depends  on  a diffuse  current.  One  model 
is  the  Hill's  vortex,  used  by  Mi  ley. 3 This  is  a rough  approxi- 
mation to  the  field  reversed  mirror.  It  is  not  self-consistent 
but  does  lend  itself  to  analytic  solutions  of  the  particle 
trajectories. 

2)  Investigation  of  electron  return  currents,  particularly  looking  al 
effects  due  to  charge  separation  after  injection  (as  suggested  by 
W.  Kunkel). 

3)  Applicaton  of  stability  codes  similar  to  that  of  Alex  Friedman's 
3~d  hybrid  code,4  using  superparticles,  to  a field  reversed  mirror 
conf i gurati on . 

b)  Study  of  transport  properties  (this  appears  to  be  a wide  open  field). 

G.  PARTICLE  ACCELERATION,  GAMMA-RAY  EMISSION,  AND  SPARKING  IN  PULSARS 

(agc  next  paqc  thfrn 

Dr.  Fawley  was  a guest  in  our  group.  This  section  consists  of 

excerpts  from  his  recent  doctoral  dissertation  (U.C.). 

3M.  Y.  Wang,  G.  H.  Miley,  "Particle  Orbits  in  Field  Reversed  Mirrors", 

CO  0-2218-508,  June  6,  1977  (submitted  to  Nuclear  Fusion). 

4A.  Friedman,  op.  ciX. 
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ABSTRACT 

A rigorous  anil  self-consistent  three  dimensional  method  is  developed  for  the  calculation 
of  steady-stale  acceleration  ol  charged  particles  in  the  open  magnetosphere  of  a rotating,  mag- 
netized, neutron  star.  This  method,  which  utilizes  the  applicability  of  electrostatic  potential 
theory  in  the  corotating  fr„me  well  within  the  light  cylinder,  succeeds  in  specifying  a unique 
current  from  pulsar  polar  caps  when  this  current  is  space-charge  limited.  When  lieid  line  curva- 
ture is  neglected,  the  only  electrostatic  acceleration  that  occurs  is  due  to  the  finite  inertia  of  the 
beam  particles.  When  curvature  is  included,  three  different  situations  may  occur  according  to 
the  geometry  of  the  system.  For  the  aligned  rotator  the  effects  of  field  line  curvature 

prevent  time-steady,  totally  non-neutral  beam  acceleration  from  the  polar  cap  surface,  whereas 
for  the  orthogonal  rotator  (/x  _ 11 ) , a time-steady  solution  does  exist  with  charged  particles  of 
both  signs  being  accelerated  from  different  halves  of  the  polar  cap.  When  the  magnetic  and 
rotation  axes  are  oblique,  half  of  the  polar  cap  surface  acts  like  the  aligned  rotator  and  half  like 
the  orthogonal  rotator.  Above  the  aligned  half,  there  is  an  electrostatically  trapped  plasma, 
while  above  the  orthogonal  half,  a space-charge  limited  beam  current  is  accelerated  to  large 
energies. 

Various  electron-positron  pair  production  mechanisms  are  studied  in  §111  and  it  is  con- 
cluded that  only  magnetic  photoabsorption  of  beam  curvature  photons  is  likely  to  create  a dense 
pair  plasma  in  the  open  magnetosphere  (unless  the  stellar  surface  temperature  is  significantly 
greater  than  100  eV  in  which  case  ion  coulomb  field  photoabsorption  or  thermal  photons  may 
produce  a large  number  of  pairs).  It  is  improbable  that  the  inertial  acceleration  corresponding 
to  the  strictly  aligned  rotator  can  lead  to  production  of  a pair  plasma  sufficiently  dense  to  break 
the  restrictive  condition  of  totally  non-neutral  current  How.  Oblique  rotators  whose  periods  are 
shorter  than  0.15  sec  /? R/ 1 7 (0,/0j):'‘ 11  cos6,  ,70,  will  have  copious  pair  production  along  favor- 
ably curved  field  lines.  Here  B | is  the  surface  magnetic  field  (assumed  dipolar)  in  units  of 
4.4 14  - 1012  gauss,  9t  is  the  actual  polar  cap  opening  angle,  0j  is  that  expected  for  a pulsar 
whose  closed  magnetosphere  extends  to  the  light  cylinder,  and  0,  is  the  latitude  of  /x  relative  to 
(1. 

In  §IV,  numerical  calculations  are  presented  for  the  x-  and  y— ray  emission  expected  from 
space-charge  limited  pulsar  models,  both  with  and  without  pair  creation-limited  potential  drops 
The  results  suggest  that  some  of  the  COS-B  localized  y-ray  scjrces  may  be  "pairlcss"  pulsars. 

A curvature-induced  acceleration  model  can  explain  the  energetics  of  the  Vela  pulsar's  y-ray 
emission  if  0jOj=2  and  5=7.  • 10i:  gauss  but  not  the  double  pulse  shape  or  the  = 140° 
separation.  For  the  Crab,  however,  it  is  not  possible  to  explain  the  energetics,  the  144”  main 
pulsc-interpulse  separation,  or  the  optical  through  y-ray  emission  spectrum  by  a surface 
acceleration  model.  It  thus  appears  likely  that  significant  particle  acceleration  must  occur  at  large 
altitudes  (Js  107  cm)  in  the  Crab  pulsar. 

In  §V,  an  analytical  solution  is  given  for  the  initial  growth  in  time  and  space  of  an  'See 
electron-positron  spark  discharge.  Such  discharges  may  occur  above  ion-cmitung  polar  cap 

v i ous 

regions  and  within  "outer  gaps".  These  discharges  can  be  sources  of  coherent  low  frequency  p.njo . 
lO(MII/.)l  radiation  in  quantities  as  large  as  10  * ergs  per  spark  in  =10  ' seconds  (Crab  pulsar 
parameters).  Finally,  in  §V|,  the  initial  results  of  an  electrostatic  simulation  of  a growing  spark 
is  presented.  These  numerical  results  confirm  the  analytical  hypothesis  that  a spark  develops 
two  ’flame  fronts"  alter  saturation  of  the  background  electric  field.  Between  these  fronts,  which 
arc  moving  rclativislically,  is  a dense  pair  plasma  with  K H=0. 
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VI.  Numerical  Simulations  of  Spark  (>rowtli  ami  Saturation 

A.  Introduction 

In  the  previous  chapter  we  gave  detailed  analytical  results  concerning  the  predicted  growth 
and  eventual  saturation  ot'  electron-positron  discharges  in  pulsar  open  magnetospheres.  We 
lound  that  electrostatic  and  electromagnetic  electric  fields  developed  in  opposition  to  the  hack- 
' ground  "vacuum"  field.  These  induced  fields  eventually  prevent  spark  growth  in  the  center  and 

we  predicted  that  the  edges  of  the  spark  would  become  relativistic  “name  fronts"  burning  their 
way  through  the  electrostatic  energy  of  the  open  magnetosphere. 

After  useful  discussion  with  FT  Scharlemann,  J.  Arons  and  C.K.  Birdsall,  the  author 
decided  to  construct  a numerical  simulation  of  sparking  to  determine  whether  the  above 
scenario  had  any  link  to  reality.  There  immed'ately  came  a choice  as  to  whether  to  attempt  a 
full  electromagnetic  treatment  or  a far  simpler  electrostatic  one.  After  study  of  the  I Ml  code 
in  use  at  Berkeley  and  Lawerence  Livermore  Laboratory.  I decided  that  it  would  be  dillicult  and 
extremely  time  consuming  to  adapt  EMI  to  handle  time-dependent  electric  fields  that  occurred 
parallel  to  the  current  J (as  is  the  situation  in  sparkin'*).  I therefore  turned  to  the  electrostatic 
code  ESI,  developed  by  A.  B.  Langdon  and  C.K.  Birdsall.  In  the  last  section  of  this  chapter  we 
discuss  the  limitations  of  a purely  electrostatic  treatment.  We  do  note  here,  though,  that  the 
* results  of  the  previous  chapter,  especially  equation  (5.45),  indicate  that  the  electrostatic 

"Coulomb"  field  may  often  dominate  the  radiation  electric  field  in  sparking  phenomena 

We  first  discuss  that  nature  of  ESI,  the  changes  made  to  it  to  study  sparks,  the  initial 

I results,  and  make  final  remarks  concerning  future  projects  for  this  type  of  simulation. 

B.  ESI  and  its  Adaptation  to  Pulsar  Sparking 

ESI  in  its  normal  form  is  a fast,  one-dimensional  (actually  l1/:  D because  it  handles  both 
v(  and  vy),  non-relativistic,  electrostatic  code  written  in  Fortran.  The  code  alternately  numeri- 
cally integrates  the  equation  of  motion  and  solves  Poisson's  equation  by  fast  Fourier  transform 
methods  for  <1*  and  E.  This  leapfrog  method  in  aceleration  — velocity  — displacement  has  been 
shown  to  be  stable  for  this  system  of  differential  equations.  The  particle  coordinates  in  phase- 
space  (x.px,p,)  are  known  "exactly"  while  the  electric  field  and  potential  are  determined  on  a 

I’  grid  system  and  interpolated  between  grid  points. 

Once  initial  conditions  have  been  specified,  the  code  simply  integrates  the  equations  forward  in 
time.  The  user  may  specify  through  input  parameters  which  diagnostics  are  to  be  printed  out 
(e.g.  phase  space  densities,  tj(.v),£(.v).  etc  ).  The  one-dimensionality  of  the  code  and  the  use 
of  cloud-in-cell  (linear  weighting  results  in  each  charge  being  "slab-like"  as  far  as  the  computer 
is  concerned.  The  slab  shape  in  .v  is  triangular  with  width  2A.v,  where  A.v  is  the  distance 
between  grid  points. 

i)  Electric  Field  Solution 

As  has  been  emphasized  throughout  this  thesis,  pulsar  electrodynamics  arc  a multi- 
dimensional problem  Thus,  it  might  seem  strange  that  we  plan  to  work  with  a I -spatial  dimen 
sion  code.  However,  it  is  important  to  realize  that  one  can  simulate  asisymmetnc  systems,  such 
as  infinite  cylinder  geometries,  by  changing  Poisson’s  equation  from 
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with  «3  = l/w 3 where  a is  the  effective  trans\erse  (i.c.  radial)  scale  of  the  system.  Here  x is  the 
coordinate  along  the  infinite  dimensions  of  the  system.  The  justification  for  expression  (6  Ih) 
follows  from  the  use  of  Fourier  transform  convolution  relationships  and  Green's  I heorem. 

The  Green  function  for  an  infinite  cylinder  of  radius  <7j,  whose  conducting  walls  are  at 
zero  potential  is 
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where  <L,  is  the  polar  cap  radius  and  jml  are  the  zeroes  of  the  Bessel  function  Jmt.  Then,  ihc 
potential  4>(.v)  results  from  convolving  the  Green  function  with  the  effective  charge  density 
i)(r')-i7*(r  ).  In  Fourier  space,  we  therefore  have 

47(A)  = £ G,  f 77(A)  — 77* (A)  ] 

(to  obtain  77(A)  and  r)*( A)  we  assume  that  they  both  go  to  zero  at  ±°°  and  are  l.e.xbesguo 
integrable).  We  now  assume  both  strict  axisymmetry  and  that  77 ( r) -r)*(r  ) does  not  vary  in 
the  transverse  &»  direction  Ue.  each  charge  is  disk-like).  We  also  approximate  the  infinite 
series  by  a single  term,  exp  — |x— .v'|/w1  , whose  Fourier  transform  is  (A3+k3)  Thus, 

(A3  + k3) 4>( A ) - +477[t,(A)-t,*(A)1  . 

which  is  exactly  equivalent  to  expression  (6  1b).  In  our  adaptation  of  HS1.  k:  was  set  to  w,  3. 
The  effects  of  this  two  dimensionality  (and  the  conducting  walls)  is  to  make  all  electric  fields  in 
the  x -direction  due  to  a point  charge  decay  in  a with  an  exponential  scale  height  of  <u, 

ESI  nominally  works  with  periodic  boundary  conditions  of  4>(0)  =<!>(/.),  where  /.  is  the 
length  of  the  system.  While  in  principle  it  is  possible  to  change  these  boundary  conditions  to 
represent  the  true  situation  of  ‘l» ( A ) — 4>(0)^0,  the  author  decided  it  would  be  computationally 
expensive  to  do  so.  Instead,  we  made  the  system  sufficiently  long  that  over  the  entire  simula- 
tion time,  no  particles  would  get  within  a couple  (L  of  the  edge  where  the  boundary  conditions 
will  affect  the  computed  induced  £-ftcld.  With  this  safety  margin,  the  induced  electric  fields  due 
to  charge  separation  within  the  spark  will  be  accurately  computed.  Most  runs  were  done  with 
16<o,  (see  § C and  Figures  6-1  and  6-2  for  additional  justification  of  this  scheme).  Rather 
than  introduce  a background  of  immobile  "pseudo"  electrons  or  positrons  to  represent  77*  and 
thus  the  vacuum  electric  field,  the  code  was  modified  to  include  a constant  to  the  calcu- 
lated induced  electric  field  This  constant  is  an  approximation  of  the  (r/R%)  -1  depen- 

dence found  above  ion-zones  and  beyond  null  surfaces.  Numerically,  was  set  equal  to 
3(1 


FIG.  6-1  Calculated  induced  Coulomb  electric  field  for  charge  density  increasing  linearly  with  a 
and  boundary  conditions  <l>(0)=*l>(  /. ).  Here  /. -Sri^  and  one  sees  that  these  boundary  condi- 
tions have  little  effect  on  the  computed  field  (in  terms  of  its  flatness)  beyond  2«,  from  the 
boundaries 
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ii)  Equation  of  Motion  and  Radiation  Reaction 

The  large  particle  energies  in  sparking  meant  that  liSI  had  to  be  changed  to  a relativism 
code.  Rather  than  keep  track  of  particle  velocities,  one  now  follows  the  particle  momenta  To 
compute  particle  movements  during  time  steps,  we  set  their  velocities  to  -*-cif  \p/nu\  '•  15  and 
compute  velocities  exactly  for  |p/wc|<15.  In  general,  only  an  extremely  small  fraction  of  the 
particles  have  momenta  below  this  limit  so  that  the  CPU  time  required  for  the  momentum  — 
velocity  conversion  is  minimal. 

The  electrons  and  positrons  follow  curved  paths  and  will  suffer  energy  losses  via  curvature 
radiation.  At  the  time  of  this  writing,  the  following  scheme  was  used  to  include  the  effects  of 
this  braking.  The  radius  of  curvature  />  is  assumed  to  be  constant  with  v and  set  equal  to 
,/3Wt . The  particle  momentum  is  integrated  forward  in  time: 

p,L  = Pou  + (6-3) 

where  p is  the  momentum  measured  in  units  of  me.  Then  two  new  quantities  arc  computed: 
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Pne«  *=  max  |/V»  ' . /’*> ) 

(6.4c) 

The  quantity  (/>,"»  - />,„.»')  represents  the  energy  lost  in  a time  interval  A/  by  a particle  coasting 
along  the  curved  magnetic  field.  The  second  quantity,  approximately  represents  a 

particle’s  radiation  reaction-limited  momentum  in  a constant  electric  field.  This  approximation 
is  quite  good  (see  Figure  6-3).  By  taking  the  maximum  of  ' and  we  insist  that  a parti- 
cle cannot  radiate  any  more  energy  than  it  would  in  reality  in  a time  interval  A/.  This  prevents  a 
particle  which  enters  a region  of  low  electric  field  from  suddenly  being  artificially  forced  to  a 
low  "radiation  reaction-limited"  energy. 
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iii)  Pair  Creation 

Sufficiently  energetic  curvature  photons  will  be  absorbed  by  the  magnetic  field  to  form 
electron-positron  pairs.  In  our  adaptation  of'  ESI,  we  assumed  a monochromatic  emission  spec- 
trum at  the  critical  photon  energy  of 

* Mic2  = SI'3  h l'i>  1 (6.5) 


and  a photon  emission  rate  of 

fl'v<  = 2 JY 
ih  9/i  he' 


4 


(6.6) 


If  the  photon  energy  computed  in  (6.5)  is  too  low  for  absorption  to  occur  in  a distance  K = 2ui,, 
then  pair  creation  was  ignored  for  the  particle  in  iiuestion.  (In  reality,  pair  creation  is  possible 
for  a distance  \~r/3',  thus  the  spark  spreads  outwards  in  v faster  than  either  our  numerical 
treatment  or  the  analytical  approach  of  the  previous  chapter  would  predict). 

An  "on-the-spot"  approximation  was  made  to  pair  creation  when  it  did  occur;  i.e.  the  pair 
was  created  at  virtually  the  same  position  as  its  parent  particle.  In  the  latest  version  of  ISI- 
SPARK,  we  take  into  account  the  "reversing"  effects  of  the  electric  field.  Thus,  one  species  is 
created  at  the  same  position  as  the  parent  particle  while  the  other  species  has  an  average  lag  of 
8x  — cAr/2.  Each  newly  created  particle  is  given  an  initial  momentum  of 


PpO'r 


Pncn  Pm' m 


cIN, 

2 — — - A/ 
dt 


(6.7) 


where  takes  into  account  the  small  pitch  angle  synchrotron  losses  of  the  newly  created 
pair 

We  also  instituted  a cutoff  to  the  pair  creation  rate  such  that  no  more  than  eight  pairs  per 
particle  could  be  created  per  time  step.  This  cutoff  is  only  needed  when  instabilities  such  as 
(OpAr— 0(1)  threaten  to  prevent  a simulation  from  finishing  due  to  exponential  blowups  One 
can  normally  tell  when  this  blowup  occurs  from  the  output  diagnostics  (/.<>.  the  pair  creation 
rate  suddenly  hits  the  maximum  of  eight  particles  per  time  step  and  the  plasma  frequency 
begins  climbing  exponentially  in  the  saturation  regime). 

ESI  normally  handles  a constant  number  of  particles.  The  growth  in  our  particle  number 
due  to  pair  creation  meant  that  another  change  had  to  be  made.  This  was  to  introduce  the  rou- 
tine "SQUEEZ”  which  squeezes  the  particle  number  down  by  an  arbitrary  factor  tin  our  case 
three  was  chosen).  Most  of  our  runs  were  done  with  a maximum  number  of  2500  electrons 
and  positrons  to  be  held  in  memory.  When  this  number  was  exceeded,  SQUEEZ  was  called  and 
it  randomly  threw  away  two-thirds  of  the  particles.  The  remaining  third  had  their  charge  and 
mass  increased  by  a factor  of  three.  Thus,  both  the  charge  density  and  plasma  frequency 
remained  unchanged,  save  for  the  random  fluctuations  due  to  the  non-continuous  distribution 
of  —800  particles  in  64-128  grid  cells.  No  problems  were  encountered  in  using  this  method  Its 
CPU  time  is  relatively  short  as  it  involves  only  one  "do-loop"  and  no  computation  ( tiist 
rcshullling  of  the  momenta  and  position  arrays). 


T 

iv)  Input  and  Output  Diagnostics 

At  the  time  of  this  writing,  FS1-SPAKK  required  the  following  input:  number  of  gud 
points(NCi),  length  in  polar  cap  radii (LTUHE),  limestep  (DT)  as  a fraction  of  d\/c,  pulsar 
period(P),  magnetic  field  in  gauss(l)FIFLD),  the  ratio  of  W,  to  WjIWCAl*.  normally  set  to  I), 
and  the  charge  (in  electron  charges)  of  the  initial  super-electron  anil  super-positron  as  a frac- 
tion of  W.T/^r/v  ai.MPFRVAL,  normally  set  to  ^MO  or  less).  The  code  then  computes  ihc 
actual  values  of  ,<!>, ,/.,  </.v,  and  dt  used  in  integration 

FS1  is  set  up  to  produce  "snapshots"  of  variables  like  charge  density  at  user-specified 
intervals  during  the  running  of  the  code.  In  our  runs,  we  found  it  most  useful  to  have  frequent 
snapshots  of  /.  and  tj,  and  less  frequent  snapshots  of  Fip)  vs  /»,  ) vs  ,v,  and  «l>.  Ai  the 

end  of  the  run,  time  histories  of  total  particle  number,  kinetic,  electrostatic,  drift,  and  thermal 
energies  are  produced. 

v)  Simplified  Flow  Chart  of  ES1-SPARK 

a)  Compute  the  induced  Coulomb  f- field  from  cloud-in-cell  weighting  of  charge  density 
using  Poisson’s  equation  modifed  for  axisymmetric  infinite  cylinder  geometry.  Add  ine  vacuum 
electric  field. 

b)  Compute  the  new  momenta  of  particles  from  equation  of  motion,  include  the  effects  ol 
radiation  reaction. 

c)  Compute  pair  production;  put  new  particles  into  momentum  and  position  array. 

d)  If  too  many  particles  to  fit  in  the  allocated  arrays,  call  SQUFEZ  and  reduce  particle 
number. 

e)  Compute  particle  velocities  and  move  particles;  compute  new  cloud-in-ccll  charge  density 

f)  Print  out  any  required  snapshots  anil  go  to  (a)  to  begin  new  cycle. 


C.  Results 

Figures  6-1  and  6-2  illustrate  the  effects  of  the  <1>(0)=<J»( A ) boundary  conditions  on  the 
induced  electric  field  from  a charge  density  linearly  increasing  with  position  in  the  \ direction. 
This  simulates  the  increase  of  tj*  with  v due  to  the  magnetic  field  curvature.  Defining  A\,„,  as 
the  distance  one  must  go  from  the  boundary  .v— 0 or  x—l.  for 
Emiiwti ( computed ) = 1 / 2 ( theoretical ) , we  see  that  A.v/Z.  drops  from  0.2  at  /.  - Soj,  to  0 064 
at  Z.=32a>l.  The  flatness  of  computed)  in  the  central  region  reflects  the  exponential  de- 

cay with  ,vof  the  <f>=0  boundary  conditions. 

Figure  6-3  plots  the  kinetic  energy  versus  time  of  a single  particle  accelerated  from  rest  by 
the  background  electric  field.  The  saturation  in  energy  occurs  because  of  the  effects  ol  radiation 
reaction.  The  computed  asymptotic  energy  of  3.73  10  compares  with  the  theoretical  radiation 
reaction-limited  energy  of  F„=3. 90-  107  [the  difference  is  due  to  our  approximations  (6.4a)  and 
(6.4b)]. 

Figure  6-4  shows  the  growth  of  a single  spark  initiated  by  a single  electron  and  positron 
pair  created  by  an  external  y-ray.  One  can  see  that  the  growth  is  initially  exponential  and  turns 
over  to  a nearly  linear  growth  rale  alter  electric  field  saturation  (see  below).  Ihc  measured  <•- 
folding  time  of  1.20  10  ’ seconds  compares  favorably  with  the  analytical  predicted  tune  ol 
2.05-10  7 for  the  chosen  parameters  iP-  0.0331 ', // =2.  • lO'-.v,  (d, /»,,)  ~ 1,  /„„  I ‘>ss/.i)s 


<> -•’  Computed  I orcnt/  factor  of  .111  electron  starling  from  rest  in  a constant  electric  field  as  a 
lunclion  ol  tune  Note  that  the  ordui ale  is  a logarithmic  scale  I lie  particle  obtains  a constant 
asymptotic  cnetgy  due  to  radiation  reaction  along  curved  tieid  lines. 
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[machine  units],  1 . 1 96/206  [statvolts/cm  |.  The  measured  rate  is  somewhat  faster  because  the 
code  permits  pair  creation  by  particles  before  they  have  obtained  radiation  reaction-limited 
energies  whereas  the  analytical  results  included  a forced  lag  time  of  r = F „/(<72/ me)  [<  /.  equa- 
tion (5.12)]. 

We  now  turn  to  the  most  important  results  of  these  numerical  simulations:  what  happen-, 
in  the  saturation  regime  of  spark  growth??  f igures  6-5  through  6-1)  presents  various  snapshots 
of  the  charge  density  and  induced  Coulomb  electrostatic  electric  field  caused  by  the  relative  dis- 
placements of  the  electron  and  positron  densities.  One  observes  that  the  electric  field  is  grow- 
ing exponentially  in  lime  and,  at  the  center  of  the  spark,  is  opposed  in  sign  to  that  of  the  back- 
ground vacuum  electric  field.  As  one  leaves  the  center  and  approaches  the  system  boundaries, 

,.,/  changes  sign.  This  occurs  us  one  crosses  the  region  of  maximum  charge  density  which 
acts  like  a surface  charge  although  its  thickness  is  comparable  to  w before  saturation,  liventu- 
ally,  Eindmnl (x=L/2)  = and  the  exponential  growth  of  the  spark  saturates.  This  saturated 

electric  field  maintains  a constant  value  and  "diffuses"  outwards  from  the  center  with  increasing 
time.  The  charge  density  docs  not  approach  t;*  but  instead  has  dr)/ dx~dr\Kl  dx  which  produces 
the  same  elTect  of  This  should  be  expected  because  the  net  charge  density 

must  be  zero  in  the  spark  and  all  quantities  will  be  symmetric  about  x=L/2. 

As  time  increases,  the  center  region  remains  shorted  out  while  the  outer  regions  have  the 
"flame  front"  structure  hypothesized  in  the  last  chapter.  These  fronts  are  topologically  similar 
to  shock  waves,  /.<».,  if  one  transforms  to  the  frame  comoving  with  a front,  all  quantities  appear 

steady  in  time.  Nearly  all  the  new  pairs  are  created  in  the  flame  front  regions  as  shown  in  I ig 

ure  6-10.  Since  the  plasma  inside  these  fronts  still  has  quite  high  Lorentz  factors  [O(10")|,  a 
small  amount  of  pair  creation  continues  to  occur  in  the  central  region  Ignoring  the  numerical 
difficulties  encountered  with  an  ever-growing  number  density,  the  code  could  run  forever  until 
the  fronts  hit  the  regions  affected  by  the  boundary  conditions. 

The  growing  number  density  and  the  slowly  "cooling"  (due  to  curvature  radiation  losses) 
pair  plasma  lead  to  the  plasma  frequency  increasing  with  time.  A relativistic  hot  plasma  has  an 

effective  to;,oc  n (AD*1  where  AT  is  the  dispersion  in  the  pair  plasma's  Lorentz  factors  Thus, 

as  n increases  and  AT  decreases,  w^Arcan  approach  I which  presents  problems.  A non-plasma 
physicist  using  this  type  of  numerical  code  quickly  learns  to  his/her  dismay  that  gross  numerical 
instabilities  ensue  because  the  imaginary  part  of  the  effective  numerical  plasma  frequency  is 
proportional  to  the  cube  of  co,,Ar.  In  Figure  6-1^.  we  show  the  beginning  of  this  instability  as 
shown  by  its  effects  upon  the  induced  Coulomb  field.  To  date,  the  only  sure-fire  method  to 
postpone  this  instability  while  retaining  the  physics  is  to  decrease  the  time  step  At.  On  occa- 
sion, the  plasma  responds  to  this  instability  by  heating  itself  up.  thus  increasing  Al',  reducing 
cjp,  and  making  the  instability  go  away  for  a few  more  time  steps.  Ikil  in  the  end.  the  instability 
wins  out,  particles  arc  accelerated  to  high  energies,  more  pairs  are  created,  n increases  faster 
than  (AT)  1 decreases,  and  there  is  utter,  complete  chaos;  the  C iiiierdiimmenun;  of  a numerical 
spark. 
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FIG.  6-4  Initial  exponential  growth  in  total  panicle  number  followed  by  saturated  linear  growth  ol 
an  electron-positron  spark  discharge. 
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CHARGE  DENSITY  at  TIME  =1-16  ns 


FIG.  6-5  Computed  charge  density  (in  machine  units)  as  a function  of  v after  I I6»iv  Tins  is  still 
in  the  exponential  growth  regime  of  the  spark.  Simulation  parameters  are  /’-  0.033 1 sec. 
B~2.  ■ 1012  gauss.  9j8a-\,  PER\'AL  = \ A0  \ fl.=  IObcm. 
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ELECTRIC  FIELD  at  TIME  =1-16  /xs 


6-6,6-7,6-9  Computed  electric  field  (in  machine  units)  due  to  spark  charge  density  separa 
tion.  Mere  the  spark  has  entered  the  saturated  growth  regime.  The  background  electric  field  i- 
+2.45  • 10'.  Thus  the  sum  of  the  induced  Coulomb  field  and  the  background  field  is  zero. 
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ELECTRIC  FIELD  at  TIME  =2-32  /xs 


FIG.  6-7 
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FIG.  6-9 
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NEWLY  CREATED  PAIRS  at  TlME=3-48/is 


6-10  The  number  density  (in  arbitrary  units)  of  newly  created  pairs  (in  the  last  time  step) 
plotted  against  .v.  Comparison  of  this  figure  with  Figure  6-9  shows  that  nearly  all  the  pairs  are 
creaif-d  at  the  edge  of  the  spark  where  K B appreciably  departs  from  zero. 
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0-  Future  Projects  in  Simulation 

I 

There  are  several  additional  physical  situations  that  should  he  simulated  besides  •>• 
infinite  cylinder  cases. 

The  original  Kuderman  Sutherland!  1975)  model  of  a higi.  ion  work  lunui  >n  polar  e.ii 
sans  space-charge  limited  How  can  be  easily  simulated  Mere  one  would  ao-lly  use  d.c  '»  o 
boundary  of  the  grid  as  a simulated  stellar  surface.  Rather  than  inserting  / , uriiia>iP>  iiu  | 

would  use  pseudo-electrons  (11  H<0  in  an  KS  model)  to  simulate  the  vacuum-induced  .*  j 

celeration.  One  must  make  the  grid  length  sufficiently  long  that  no  particles  will  get  within  i , 

few  pf  the  outer  boundary  .v  = L over  the  course  ol  the  run  One  crucial  question  imuetmnp  j 

the  original  RS  model  is  whether  a stable  "vacuum  gap”  actually  does  form  above  the  stellar 
surface,  or  is  an  outward  moving  flame  front  (behind  which  r/  relaxes  to  r)A.  all  the  way  down  to 
the  stellar  surface)  more  likely?  Given  that  cross-held  marching  inhibits  an  gap  front  obtaining 
an  equilibrium  configuration,  the  flame  Iront  picture  nta>  he  more  realistic 

i 

Another  problem  is  to  follow  the  inward  moving  flame  front  in  a null  -surface  tvpe 
configuration.  If  one  can  somehow  include  a stellar  surface  twh  induced  space-charge  limited 
flow),  then  it  may  be  possible  to  establish  what  is  the  final  equilibrium  of  an  "outer  gap"  spark 
discharge  (admittedly  ignoring  cross-field  drifts)  1 he  autnor's  present  opinion  is  that  the 
equilibrium  may  be  similar  to  that  shown  in  Figure  b-li.  There  will  be  a dense  plasma  extend- 
ing from  the  stellar  surface  to  the  surface  rni-.x<t)>.«f>,)  where  rmj,  is  the  maximum  radius  at 
which  pair  creation  will  still  occurs.  Beyond  there  is  a totally  non-neutral  space-charge  lim- 
ited current  flow  (which  may  balance  the  current  ol  the  cos  >0  hall  ol  the  open  magneto- 
sphere). T his  current  may  be  able  to  extract  a significant  fraction  at  the  potcnii.il  drop  available  ; 

in  a totally  vacuum  open  magnetosphere.  II  so,  this  type  ol  model  might  be  more  successful  in 
explaining  the  energeties  of  the  Crab  pulsar's  emission  than  are  surface  acceleration  models 
such  as  ours  or  that  of  Kuderman  and  Sutherland.  | 

Finally,  the  author  hopes  to  do  a more  accurate  electromagnetic  simulation  ol  sparks.  I his 
will  involve  making  major  adaptations  to  F.Ml.  It  may  also  require  the  introduction  ol  the  con- 
ducting boundary  surface  between  the  open  and  closed  magnetosphere.  Thus,  the  problem  may 
become  more  of  a wave-guide  one.  It  is  doubtful  that  any  of  this  will  be  easy.  Unless  the  con- 
ducting walls  introduce  peculiar  elTects,  it  would  be  surprising  if  the  electromagnetic  eflects  rad- 
ically change  the  nature  of  the  flame  front  equilibria  found  in  our  electrostatic  simulations 
Based  upon  our  analytical  results,  we  suspect  that  the  radiation  electric  field  is  only  important 
just  before  electrostatic  saturation  occurs  in  the  center  of  the  spark;  the  electrostatic  field 
saturation  might  be  postponed  lor  a few  e-folding  times  but  it  will  eventually  always  dominate 
the  radiation  field. 


FIG.  6-12  Possible  equilibrium  situation  of  an  n ft<0  oblique  rotator.  Spark  discharges  occur  in 
both  halves  of  the  open  magnetosphere.  The  outward  moving  flame  fronts  have  stopped  at  a 
region  where  pair  creation  is  barely  possible  because  of  the  fallofT  of  the  background  dipolar 
magnetic  field.  Beyond  this  region,  there  are  non  neutral  beams  of  electrons  and  positrons  in 
the  cos</>,<0  and  cos</».>0  regions  respectively. 
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Se.cti.un  171 

COOK  DliVliLOPMI-NT  AND  MAINTKNANCI; 


A.  ESI  CODE 


Nothing  to  report  this  quarter. 

B.  EMI  CODE 

Nothing  to  report  this  quarter. 

C.  ELECTROSTATIC  2*d  CODE  DEVELOPMENT  (EZOHAR) 

W.  M.  Nevins  and  Dr.  Yoshi  Matsuda 

We  have  developed  a new  version  of  EZOHAR  which  employs  a 
guiding-center  (ExB)  mover  for  electrons  and  the  usual  leap-frog  mover 
for  ions.  The  purpose  is  to  simulate  drift  instabilities  such  as  drift 
cyclotron  and  drift  cone  modes.  The  guiding  center  mover  and  its  algor- 
ithm used  here  was  described  earlier  in  QPR  III,  1977- 

A number  of  tests  runs  have  been  performed  to  check  effects 
of  boundary  conditions,  stochastic  heating,  and  frequency,  and  frequency 
spectrum.  We  have  used  a uniform  plasma  in  the  bounded  system  described 
earlier.  The  test  results  show  that  the  density  stays  uniform  during  the 
entire  simulation  run,  typically  1000  steps,  causing  no  apparent  ill 
effects  near  the  wall.  The  stochastic  heating  observed  is  shown  in 
Fig.  1.  The  total  energy  increase  in  this  case  is  about  0.3£  in  1000 
steps,  which  is  reasonably  small.  Figure  2 shows  measured  frequency 
spectra  for  cyclotron  harmonic  waves.  These  figures  are  for  two  differ- 
ent wave  numbers  corresponding  to  k p.  0.1*  and  2.8  where  i>.  is 

r J y i i 

ion  Larmor  radius.  We  can  clearly  see  the  cyclotron  harmonics  as  expected. 

These  results  indicate  that  the  guiding-center  version  of 
EZOHAR  is  working  properly.  We  will  be  using  this  code  to  simulate  drill 
cyclotron  instability  and  study  its  linear  and  nonlinear  behavior  in  the 
comi ng  quarter. 


nergy  ( ar  bitrary  unit  ) 
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D.  THE  COLD  BEAM  INSTABILITY;  PREDICTION  OF  SATURATION  LEVEL* 

W.  M.  Nevins  and  J.  Albritton  (visiting  from  Univ.  of  Rochester 
Fall,  1976) 

In  the  limit  of  a fast  cold  beam  B t v /u>  Ax  » 1 the  arid-alias 

op 

cold  beam  instability  may  be  understood  as  resulting  from  the  interaction 
between  the  beam  (whose  velocity  is  v ) and  the  p=l  alias  of  the 


plasma  wave,  which  has  a frequency 


w = k.v 
I o 


where 


* (k‘pkg) 


9'p-l 


k = 2n/Ax 

g 


In  fact,  the  growth  rate  of  this  instability  may  be  obtained  by  keeping 
only  the  p = 0 and  p=l  terms  in  the  finite-grid  dispersion  relation1 
to  obtain 


1 - w (k) 
P 


(u>  - kv  ) ‘ 
o 


(w‘klVor 


We  have  used  linear  weighting,  a three  point  Laplacian,  and  two  point 

2 

centered  differencing  for  the  gradient.  Hence,  u>p(k)  »s  given  by 


| (k)  = oo2  dif  (kAx)  dif2  (^p) 


It  is  also  convenient  to  define  a frequency 


£2  = a)  - kv 


-Work  done  Fall  1976;  saved  for  possible  journal  note,  which  has  not  yet  been 
written.  Previous  work  was  a report  in  our  QPR’s  during  1975,  under  tin- 
title  "Cold  Beam  Instability". 

1 Due  to  A.  B.  Langdon,  given  explicitly  in  C.  K.  Birdsall,  N.  Marion,  G.  Smith, 
"Cold  Beam  Nonphysical  Instabilities  and  Cures",  Proc.  7th  Conf.  on  Numerical 
Simulation  of  Plasmas,  NYU  June  2-4  1975.  See  Eq . (5),  P-  178. 
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when  w/k , « v 

I o 


Strong  coupling  between  the  alias  mode  and  the  beam  will  occur 
i .e.  , when  ft  is  small.  Hence,  we  consider  the  limit 
In  addition,  we  restrict  out  attention  to  small  values  of  k , 

The  leading  terms  in 


ft/m  « 1 

P 

such  that  k/k  « \ , and  assume  that  k v /to  » I 
9 g o p 

the  dispersion  relation  are  then 


where 


°22  > - C3 

“p  (i-O3 


('  - '*<’) 


(4) 


Hence,  ft  is  purely  imaginary  for  small,  positive  r,  . The  maximum 

growth  rate  is  determined  by  choosing  r,  to  maximize  the  absolute  va  1 1 
2 

of  ft  . We  find 

(6) 


Y = 0.2  w 

max  p 


at 


iT=  oA 

9 


(7) 


which  agrees  with  the  more  complete  (more  values  of  p)  solutions  of  Bird- 
sal  1 zt  al.2 

In  this  limit  (i.e.,  k v /w  »1),  saturation  occurs  due  to  the 

9 o p 

trapping  of  beam  particles  in  the  electric  field  of  the  alias  wave.  This 
trapping  occurs  when 


k|6l  ~ 1 


(8) 


whore  <*>  j is  the  (lagrangian)  di  splacement  of  the  beam  particles  due  to 
the  presence  of  the  alias  wave. 


2C.  K.  Birdsall  vtaJL.,  op.  ett. 
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We  relate  the  velocity  spread  of  the  beam, 
i an  displacement  by  the  approximation 


, to  the  Lagrang- 


v = 
t 


|6,|  - |Q«)| 


(9) 


In  writing  the  second  equality,  we  have  used  the  fact  that  SI  is  the 
frequency  of  the  alias  wave  in  the  reference  frame  of  the  beam.  Since 
trapping  occurs  when  k)6J=l  , it  follows  that  the  beam  spread  at 
saturation  is  given  by  (approximately) 


v 


t 


(10) 


Substituting  the  appropriate  values  for  the  fastest  growing  mode  gives 


w. 


saturat ion 


- 0.05  ^ Ax 
P 


(11) 


Finally,  we  note  that  when  the  Debye  length  of  the  beam  is  taken  to  be 
Xp  = v^iDp  , this  result  may  be  written  as 


X 


D 


Ax 


* 0.05 


(12) 


This  is  in  excellent  agreement  with  the  simulation  results  shown  in  Fig.  I 
on  p.  18,  QPR  II,  July  1,  1975.  We  also  note  that  an  examination  of  the 
phase  space  plots  at  saturation  reveals  clumping  of  the  beam  into  vortices 
with  a spacing  equal  to  the  wave  length  of  the  p=l  alias.  Hence,  we 
conclude  that  this  simple  trapping  model  accurately  describes  the  satur 
ation  of  the  grid-alias  instability  in  the  limit  of  a fast  cold  beam. 
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Section  IV 

PLASMA  SIMULATION  THXT 

No  progress  this  quarter. 

Section  V 

SUMMARY  OF  REPORTS,  TALKS,  PUBLICATIONS  IN  PAST  QUARTER 

M.  J.  Gerver,  "Radial  Normal  Mode  Calculation  of  Warm  Plasma 
Stabilization  of  the  Drift  Cone  Mode",  Phys.  Fluids  3.  PP-  M3-M6, 
March  1978. 
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